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Abstract 

A phase-field method is applied to the modeling of flow and breakup of droplets in a T-shaped 
junction in the hydrodynamic regime where capillary and viscous stresses dominate over inertial 
forces, which is characteristic of microfluidic devices. The transport equations are solved numer- 
ically in the three-dimensional geometry, and the dependence of the droplet breakup on the flow 
rates, surface tension and viscosities of the two components is investigated in detail. The model 
reproduces quite accurately the phase diagram observed in experiments performed with immiscible 
fluids. The critical capillary number for droplet breakup depends on the viscosity contrast, with a 
trend which is analogous to that observed for free isolated droplets in hyperbolic flow. 

PACS numbers: 47. 11+j, 47.60. +i,47.55Dz,47.55.Kf 
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I. INTRODUCTION 



Much of the recent interest in the behavior of fluids in small channels stems from the 
wide range of applications that microfluidic devices are finding in biology and chemistry. 
The laminar regime of the flow in these systems allows an accurate control of mass and 
momentum transport, such that chemical species can be brought together in small quan- 
tities within restricted regions of space, enabling, for example, fast mixing and high 
resolution of reaction kinetics [2|. In the case of two-phase systems, nano- and pico-liter 
droplets can be handled with great precision when the destabilizingeffects of wall impurities 
can be suppressed, producing well controlled structures ; partial wetting, in 

fact, introduces non-linear instabilities which may generate complex erratic [t| or oscilla- 
tory patterns . The variety of two-phase flow behaviors (see figure ^) is both interesting 
at a fundamental level and relevant for those applications based on the accurate tuning of 
droplet sizes [jj: micro-droplets could, for example, be used to confine chemical reactions, 
or to deliver substances in well defined amounts. Such a control is generally attainable for a 
limited range of flow conditions and with fluids having properties that are compatible with 
the surfaces of the device used. The combinations of materials, geometries and flow condi- 
tions which are being used in experimental setups are rapidly growing, and modeling can be 
extremely helpful both in providing insight and in improving the design and prototyping of 
microfluidic devices. In this respect, the main challenge of a modeling technique is to cope 
with the three-dimensionality and the strong surface stresses present both at the fluid-fluid 
and fluid-solid interfaces. 

Phase-field models have been developed for the study of the dynamics of critical phe- 
nomena, in particular in relation to nucleation and spinodal decomposition . The spatial 
variations of the phase order parameter indicate the location of interfaces separating the 
growing domains, and interface coalescence and breakup is driven by the local curvature- 
dependent solubility (Gibbs-Thomson effect): interfacial regions of high curvature have a 
higher chemical potential which drives a mass flux to the surrounding medium to reduce the 
curvature and therefore the chemical potential gradients. The effects of fluid flow are in- 
cluded by adding an advection term to the diffusion equation of the order parameter, which 
is then coupled to a momentum transport equation to describe the hydrodynamics of a mul- 



tiphase system 
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121 ] . From a computational point of view, modeling of fluid boundaries 
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as having finite thickness greatly simplifies the handling of interface coalescence or breakup 
events, and in the case of phenomena occurring near the critical point, this approach is fully 
justified by the fact that, in the mean field approximation, interfaces are indeed diffuse, with 
a characteristic length which is comparable with the typical domain size. On the other hand, 
it can be shown that phase-field physics approaches asymptotically that of two immiscible 
fluids as the interface thickness becomes small compared to the local curvature, both for 
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unbounded flow and in the presence of solid surfaces 

In this paper we use a phase-field model to study the three-dimensional flow of a binary 
mixture in constrained geometries in the hydrodynamic regime which characterizes the flow 
in microfluidic devices, i.e. with low Reynolds and Capillary numbers; we show how the 
computed droplet flow and breakup dynamics in a T-shaped junction follows closely that 
observed experimentally for water in oil mixtures, despite the fact that the interface thickness 
in the simulations is far from the sharp-interface limit. 

The paper is structured as follows: section |n] summarizes how the thermodynamic and 
transport properties of the non-homogeneous system can be derived from a free-energy 
functional; section ITTTl describes the numerical method used to solve the transport equations 
while section IIVI presents the free-flow tests which were performed to validate the imple- 
mentation of the code. Finally, in sections IVl and fVTl the results for the modeling of droplet 
flow and breakup in a T-shaped junction are reported. 



II. MODEL EQUATIONS 

A unifying feature of all phase-field models is the existence of a free-energy functional, 
which, besides determining the equilibrium properties, strongly influences the dynamics of 
the system, whose transport equations can be derived from conservation arguments Q|, or 
from generalized hydrodynamic theory . In this section we briefly review this theoretical 
framework for a two-phase system, illustrating in particular how the surface tension, the 
wetting properties and the pressure tensor follow consistently from the choice of the free- 
energy functional. 



3 



FIG. 1: Left: examples of simulations of two-phase flow dynamics in microfluidic devices: a) 
droplet breakup in a symmetric T-junction; b) and c), droplet formation in a T-junction with 
cross flow; d) droplet formation in a flow-focusing cross-junction. Right: frames showing the 
modeled three-dimensional flow of a droplet in a symmetric T-junction. 

A. Thermodynamics of the non homogeneous system 

We consider a mixture of two fluids, A and B, with the Cahn-Hilliard-van der Waals 
form for the free energy in the volume V, 



where n = ha + ub is the total particle number density, (p = ns/n is the molar fraction of 
component B, k determines the surface tension (see section Hi Bl below) and / is the sum of 
the free-energy densities of the pure components, which for simplicity is assumed to depend 
only on the total density. Also, W is the free energy of mixing, for which we shall take the 
expression for a simple symmetric mixture 



which is the sum of the ideal and the excess parts, k is Boltzmann's constant, and T and 
T c are, respectively, the absolute and critical temperatures. The critical point corresponds 
to the molar fraction ip c = 1/2 and temperature T c , so that below T c the system phase 
separates; at T = the two components completely demix so that the equilibrium bulk 
molar fractions become <pa = and cpg = 1. In the incompressible limit, n =const., the 




(1) 



W(<p) = 2kT c ip(l -ip)+kT [(^log^ + {l-<p) log(l - tp)] 



(2) 



equilibrium concentration satisfies the Euler-Lagrange equation 

M (r) = -«VV(r) + W'(<p(r)) = /w, (3) 

where /i(r) = 5F/5(f(r) is the chemical potential difference, and /i coe x is the Lagrange multi- 
plier which fixes the total molar fraction. For a flat interface /i coe x = and (p is a function of 
one spatial coordinate only: the concentration profile varies between the equilibrium molar 
fractions ipA an d <Pb, within a region having a finite characteristic thickness £ proportional 
to a/ n/kT c . In non-equilibrium conditions, local imbalances of the chemical potential n(r) 
will generate diffusion currents which tend to restore the configuration satisfying equation 



B. Surface tension 



The mechanical definition of surface tension for a flat interface is given by the integral of 
the difference between the normal and transverse pressures: 

o = j°° dx (p„ - 5^±^ , (4 ) 

where x is the direction normal to the interface. The pressure tensor must satisfy the 
Gibbs-Duhem equation (we consider only isothermal transformations) 

V,3 • P a /3 = n A V a ^A + n B V a fi B , (5) 

where = SF/Sua and /x& = SF/Sub are the chemical potentials computed as functional 
derivatives of the free energy (JTJ). The integrand in the functional (JTJ is a Lagrangian density, 
and is invariant under space translations; Noethers theorem can be applied to construct the 



expression for the associated conserved current 



121 ]. leading to the explicit form 



P a/ 3 = 8 a p [nf{n) - f(n)) + nnVapVp^, (6) 

where we recognize the non-isotropic part (the second term in (JHJ)) first introduced by 
Korteweg. From equations (@J) and (0), we get 



DC 



a = I dx Kn [ — ) , (7) 
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which corresponds to the excess energy with respect to the homogeneous system, also called 
the distortion energy For n =constant, the integral in (JJJ) is easily expressed in terms 
of the mixing energy W: 



a = n dip \J2kW((p), (8) 
J ip A 

with W(ip) = W(<p) — W(<pa) > 0. Without loss of generality, we can put k = kT c ^ 2 , so 
that the surface tension is a = kT c n^Q, with Q given by the integral 



•rD 



fi = / dip \J2W(p)/kT c . (9) 

Jip A 

In the particular case T = 0, the one-dimensional equilibrium molar fraction is <p(x) = 
1/2 + 1/2 sin(2rr/£) for |x| < vr^/2, or 1 otherwise. The interface width is 7r£ and the surface 
tension integral Q gives Q = it/ A. For T > the bulk concentrations p>A and p>B of the flat 
interface profile are approached asymptotically, and the exact definition of £ can be assigned 
in an arbitrary way. In this work £ is defined as the interval over which the flat interface 
concentration rises from 34.4% to 65.6% of the total variation cpg — p>A > 0. In other terms, 
given the profile ip(x), £ = \x2 — xi\, where (<p(xi) — <pa)/(<Pb — <Pa) = 1/2 — 1/2 sin(l/7r), 
{<p(x2) — Pa)/{.Pb — Pa) = 1/2 + 1/2 sin(l/7r). In the simulations reported here we chose 
T/T c = 0.8, such that the mixing energy (J2J) is far from the singular limit T = 0, and the 
interface thickness is still quite small. 



C. Wetting 

In the framework of the Cahn-Hilliard theory of diffuse interfaces, the interaction of the 
fluid components with a wall is introduced by adding to the functional ((TJ) a surface energy 
term of the form 

F s [n A ,n B } = [ dr5(S)$(n A ,n B ), (10) 
Jv 

where 5 is Dirac delta function, S(r) = defines the locus of the points of the surface 

n 

and $ is the surface energy density. From the series expansion of $ [18|], we keep the 
first-order term only, $ = — k(Gatia + Gb«b), with the coefficients Ga, Gb weighting the 
fluid-wall energy with respect to the same parameter k governing the fluid-fluid surface 
tension of equation (JJJ). In the incompressible case, surface energy variations depend only 
on the molar fraction ip and the difference G = Gb — Ga- The addition of the surface 



6 



free energy (|1U|) modifies the Euler-Lagrange equation Q in the proximity of the wall and 
imposes a discontinuity in the gradient of (p: 

- KVV + W\ip) - KG?«J(5) = /icoex! (11) 

the wall interaction is therefore included in the model as a boundary condition by simply 
imposing the constraint 

n-Vip = G, (12) 

where n is the normal to the wall surface, pointing inward to the solid. The parameter G, 
which has dimensions of 1/Length, determines the wetting of the solid by the two phases; 
G > implies an attraction for component B, G < a repulsion. As in the unbounded case 
of equation @, the properties of one-dimensional solutions of equation (fTTj) can be deduced 
from the first integral 

-W((p) = 0. (13) 



K 
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dx 

In particular, the equilibrium values (p s of the molar fraction at the surface correspond to 
the solutions of 

\G\ = J-W{<p), (14) 

V K 

and the solid-fluid surface energies a a and <jb for the two phases are obtained as the sum 
of —nnGips and the distortion energy integrals of the from (JHj) Q, 0|- Together with the 
fluid-fluid surface tension a, the surface energies a a and <jb determine the static contact 
angle 6 through Young's law, cos 9 = (a a — ctb)/^- 



D. Transport equations 

The transport current for the density % (ha = n — Ub, since we keep n constant) is the 
sum of an advection and a dissipative term, there the latter is proportional to the mobility 
A and the chemical potential gradient V/i. The fluid momentum equation, besides the non- 
linear advection term, contains both reactive and dissipative forces, depending respectively 
on the pressure tensor and the shear viscosity 77: 

dn B 



()f V • (n B v) + V • (AV/i) (15) 

^ = -V • P - V • (pvv) + V ■ t?[Vv + (Vvf], (16) 



where p is the mass density and v is the mass average velocity. The isotropic part of the 
pressure stress tensor (|5jl. p = nf'(n) — f(n), is solved for by imposing the constraint Vv = 
rather than being derived by a specific form of the bulk free energy of the system. We will 
also take equal molar masses and a constant mobility for the two components, A = Dn/kT, 
with D the diffusion coefficient. The viscosity changes with the molar fraction, p = rj*h(ip), 
and the following functional form is assumed for h: 

f tanh (^) + 4±± if A > 1, 

h(<p) = I V 2A ' ) X ,( (17) 
\ (*=i) tanh (^) + A±l if A < 1 , 

where A = Pb/va is the viscosity ratio and the parameter x determines the sharpness of the 
transition of the viscosity as the phase separating interface is crossed; when \ ^ — <Pa, 
we have that p* = max{pA, ps} is the viscosity of the more viscous phase. 

Given the unit time T, length L, mass M, and density h, we have the following dimen- 
sionless groups: Peclet number Pe = L 2 /DT, Reynolds number Re = pVL/p* (V = L/T), 
Weber number We = pV 2 L/a and the Cahn number C = We will maintain the sym- 
bols v and p of the dimensionless fields (p is rescaled with M/LT 2 ), W for the rescaled 
mixing energy WT 2 /ML 2 , t and V for the dimensionless time and differential operator. 
The dimensionless transport equations become 

^ = -V • (^v) + i-V 2 [-C 2 W + W'{<p)] (18) 

_ = _V.(vv)-Vp-wSV.(V^) < 19 > 

+ ^-V-/t[Vv + (Vv) T ]. (20) 
Re 

The phase-field Navier-Stokes equations converge to the classical sharp interface behavior as 
the interface thickness, indicated by C, is reduced to zero along with the diffusivity 1/Pe jl^ . 
The mesh spacing used in the numerical solution imposes the main constraint on the Cahn 
number, since a sufficient number of points is needed to describe the interface profile in order 
to avoid spurious effects and to control the convergence of the solver. At finite C straining 
flows can thin or thicken the interface, and this must be resisted by diffusion, which therefore 
requires 1/Pe to be high enough. On the other hand, too large a diffusion will overly damp 
the flow, which imposes a lower limit on the Peclet number. In our simulations, C and Pe, 
rather than being referred to the realistic values of an immiscible mixture such as water in 
oil, will therefore be chosen to optimize the approximation of the sharp interface physics, 
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compatible with the limits imposed by the mesh size. Despite the fact that C will be typically 
as large as 1/20 and Pe of the order 10, which is determined by the fact that we have to 
deal with three-dimensional simulations, we will see that the model is capable of capturing 
the main features of droplet flow in simple microfluidic devices. 



III. THE NUMERICAL ALGORITHM 



The differential equations ()18|) and (j2*Uj) are discretized on a uniform three-dimensional 
Cartesian grid with staggered velocities; the molar fraction and pressure fields are collocated 
at the centers of the cubic cells, while the velocity components are placed on the faces, and 
the boundaries of the simulated domain always coincide with a face of a grid cell, be it a wall, 
an inlet or a pressure outlet. The time evolution is implemented with an fully implicit Euler 
scheme, which is first order in time. Given the array for the discretized fields at time step 
n, 0" =o 3 = ip n , v ™, v%, t>3 and p n , their values at time step n + 1 satisfy a set of equations 
of the form 

<^ +1 = <T- ^tN i {r 1 + \4> n 2 + \<t>T l AT\p n+l ), (2i) 

with the constraint of incompressibility, V • v n+1 = 0. The explicit form of the nonlinear 
functions Ns is 



N o m +1 },P n+1 )= (22) 

V ' a { V n+1 v n a +l ) - -j- [-C 2 VV +1 + V 2 W'( V n+1 )] , 
re L J 

iV Q ({C +1 },p n+1 ) = (23) 

v,« + V 1 ) + v aP n+1 + ^v,(v a ^ +1 v^ +1 ) 

-^V^MV^ 1 + V^r 1 )] a = 1, 2, 3, 

and the operators V a refer to the finite-difference approximations of the partial derivative 
operators. The iterative scheme for the completion of the time step goes as follows: 

1. Given the current guess of the new values for the fields, <fi* and p*, which at the 
beginning of the iteration loop will be simply 0" and p n , the equations (|2~T]) are solved 
in sequence to get the new guesses 0** 

cp** = <p n -ddN^<p**,vl,vl,vt,p% (24) 



vT = v?-AtN 2 &**,v?,v*,v*,p*), (25) 
v** = v2-ArN 3 (<p**,v?,v?,vlp*), (26) 
v? = v%-ATN 4 (<p*\vr,v**,v?,p*). (27) 

Each inner iteration is solved recursively by substituting 0** ~ 4>k+i = <Pk + 8k with 
4>o — n an d linearizing in 5k- 

2. the continuity equation, is enforced with a prediction-correction method (see, for 
e^ple, ref . H), leading t o new values v*** and p**, 

3. the outer iteration is completed by setting ip** — > tp*, v*** — >■ u*, p** — > p*, and 
computing the of residuals 

Ri(ti) = <l>i-<l>*i-Ni({$},P*), i = 0,l,2,3. (28) 

together with the divergence field D = V • u*; the outer iteration loop is completed 
if H-Ro^*)!! < e\\(p*\\ and ||i?j(^*)|| < e\\v*\\, i = 1,2,3, and \\D\\ < e, all satisfied. In 
all simulations presented in this paper we set e = 10 -5 . When the desired accuracy is 
reached, the new fields are obtained by setting = </>*, p n+1 = p*. 

Besides the periodic boundaries used for the case of unbounded flow, we can handle walls, 
inlet and outlets. For all three kind of boundaries the normal gradient of the chemical 
potential is set to zero (n ■ V/i = 0). The boundary conditions (BC) at the walls include 
the constraint (J12|) for the molar fraction, the no-flow conditions n • Vp = and vi = 0, 
and the no-slip condition vy = v wa n, where v wa n is the tangential velocity of the wall. Inlets 
have Dirichlet BC for the molar fraction and velocity and Neumann BC for the pressure 
{ip = const., vj_ = const., vy = 0, n • Vp = 0), while outlets have the complementary setup 
(n ■ V<£> = 0, n ■ Vv = 0, vy = 0, p — const.). 

IV. FREE FLOW DYNAMICS 
A. Spurious currents 

The appearance of spurious currents in the proximity of the phase boundaries is a common 
feature shared by all numerical methods which include capillary stresses in the Navier- 



Stokes equation as a volume force to be computed on a grid 
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22j . A simple numerical test to 



FIG. 2: Visualization of spurious currents for a static droplet relaxed on three different cubic 
meshes: a) Ax = C, b) Ax = 2/3C, c) Ax = 1/2C, C = 0.1. The pressure, concentration and 
velocity fields have reached the stationary state, and parasite currents persist in the proximity of 
the interface region indicated by the concentric circles. The spurious currents reflect the four-fold 
symmetry of the lattice, and depend on the size of the mesh. Measured in units of the reference 
velocity Re/We, the maximum spurious velocity decreases from 4 • 10 -3 (a) to 2 • 10~ 3 (b), and 
3 • 10~ 4 (c). 

illustrate their characteristics is to consider a droplet at rest; when the pressure and molar 
fraction fields are equilibrated, the velocity field is not zero, and parasite currents appear 
in the interface region, with a pattern which reproduces the symmetry of the underlying 
lattice (figure |2J). The origin of the spurious currents can be explained by analyzing how 
the discretized form of the differential operators modifies the continuum transport equations 
(see Appendix), and their magnitude grows with the coarseness of grid and the ratio between 
surface tension and viscosity. In our simulations the size of spurious currents was controlled 
by changing the lattice spacing, keeping the modulus below 5% of the characteristic velocity 
of the flow. 

B. Droplet deformation and breakup 

Droplet deformation and breakup in unbounded simple shear flow depends on the ratio 
between droplet and matrix viscosities, A = rja/rjc, and the capillary number Ca = 7] c jR/a 
where 7 is the shear rate and R is the radius. For low Ca, surface stresses are balanced by 
viscous ones and the droplet is stretched to a stationary ellipsoidal shape, with its principal 
axis following the direction of maximum elongation; the deformation is expressed in terms 
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of Taylor's parameter D = {L — B)/(L + B) 23|, L being the long breadth of the drop, 
and B the small breadth (see inset to figure EK)- The angle 9 between the axis of maximum 
elongation and the shear direction gives the orientation of the droplet. Simulation results 
are consistent with small deformation theories showing a linear relation between D and 
8 and Ca for Ca <C 1 (figure Eb). As the capillary number is increased above a critical 
value Ca c , no steady shape exists and the droplet breaks. In our simulations, we find 
0.4 < Ca c < 0.44, in agreement with the value Ca c ~ 0.41 obtained with more accurate 
numerical methods 124, 251. 





FIG. 3: Variation of deformation (a) and tilt angle (b) as a function of capillary number for a 
droplet in simple shear. The straight lines show the Ca <C 1 limits derived by Taylor 23] and 
Chaffey and Brenner [2^ (A = 1). The simulations were performed on a 80 x 60 x 60 box with the 
top and bottom faces moving with opposite velocities parallel to the longer side. C = Ax = 0.05. 



V. DROPLET FLOW IN SQUARE MICROCHANNELS 

As a first example of two-phase flow in a constrained geometry we discuss the the pres- 
sure driven motion of droplets inside a square channel. By choosing the appropriate value 
of the parameter G (equation (|TTj) h the contact angle to the wall is set to zero, so that the 
continuous phase completely wets the boundaries. Exact results for droplet motion in cap- 



illaries are known only in the case of circular tubes 



In this case, theory shows that the 



speed of the drop U exceeds the average fluid speed by and amount UW, where W oc Ca^ 3 , 
and CaT = VcU/cr is the capillary number. In the case of square channels, to our knowledge, 
theoretical and experimental studies have been published only on the flow of bubbles, that 
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29l I30L |3l| . For viscous drops in square 



is in the limit of small droplet viscosity r)d — > |28 . 
channels we verified that the value of A strongly influences the droplet speed; at a given Cax 
the fractional velocity correction W increases as the viscosity ratio is decreased. We tested 
in particular the motion of droplets whose diameter is comparable with the channel width 
(figure EJ), and verified that the values of W in the range 0.01 < Cax < 0.2 are consistent 
with the experimental results for bubbles by Thulasidas et al. i.e. 0.3 < W < 0.45. 
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FIG. 4: Droplet flowing in a capillary with square cross section of width w; the (normalized) 
droplet excess speed is plotted against the capillary number for different viscosities: A = 1/32 (O)i 
A = 1/8 (A), and A = 1/2 (□). U is the droplet speed while the V is average flow speed of the 
continuous phase (see inset). The undeformed radius of the droplet is 0.55w. The simulation mesh 
is made by 125 x 25 x 25 points, C = Ax = 0.04. 



VI. DROPLET BREAKUP IN A T JUNCTION 

These simulations were motivated bv the experiments by Link et al. on droplet breakup 
in a T-shaped junction configuration 6]. Isolated water droplets dispersed in an oil matrix 
are transported by a pressure-driven flow along a hydrophobic microchannel until they 
reach the junction, where they are stretched by the local elongational flow. Depending on 
the droplet size and relative strengths of viscous to capillary stresses, they split or simply 
flow down one of the two branches (see insets a) and b) in figure EJ). The dimensionless 
groups are the capillary number Cax, the viscosity ratio, A = r]d/r] c and the initial droplet 
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extension e = lo/nwo, where Iq and u>o are the length and the diameter respectively of the 
droplet upstream of the T-junction. 

Figure summarizes the results for two values of A by reporting for each simulation a 
point with coordinates (eq, Cax); empty and filled symbols mark breaking and non-breaking 
droplets respectively. These two regions are separated by a curve curve whose shape can be 
derived with a general argument based on the Rayleigh-Plateau limit j^|, and is given by 
the function 

Ca c = ae (l/e 2 /3 - l) 2 , (29) 

where a(A) is a constant which depends on the viscosity contrast. We verified that this 
function describes very accurately the critical line for a range of A. In particular we get 
a(l/8) = 1.02 ± 0.03, a(l/2) = 0.60 ± 0.04, a(l) = 0.50 ± 0.01, a(3/2) = 0.35 ± 0.04, 
a(2) = 0.33 ± 0.04. These data clearly indicate that droplet breakup in the T-junction 
occurs more easily as the viscosity ration is increased, even when A > 1. The variation is 
quite rapid as the region around A = 1 is crossed, reaching a uniform value for large A. Such 
behavior is in agreement with the viscosity dependence of the critical capillary number for 



a free droplet in a purely elongational flow 



321 ] . Our result for a (1/8) closely agrees with 



the experimental value (a = 1) obtained with water droplets dispersed in hexadecane oil ^| 
(viscosity ratio A = 1/8), and is a demonstration of the accuracy of the numerical method 
presented here. As a final remark on the computational side, we would like to point out 
that the position of the critical line remains quite consistent for different values of the Cahn 
number; figure El in fact reports data obtained for C = 1/15 and C = 1/20, showing how 
such relatively large values are sufficient to reproduce the critical boundary described by the 
function (J2HJ). 



VII. CONCLUSIONS 



We have discussed the application of a phase-field method to the modeling of two-phase 
flow in microfluidic devices in three dimensions; in particular, this approach is validated by 
investigating the dynamics of droplet breakup in a symmetric T-shaped junction. The phase 
diagram obtained from the numerical simulations agrees very closely with the experimental 
results which were performed with water-in-oil dispersions, confirming the reliability of the 
numerical approach. This conclusion was not guaranteed, given that the physics of the 
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FIG. 5: Bottom: (eg, Cax) diagram summarizing the conditions for breaking drops at a symmetric 
T-junction for A = 1/8 (left) and A = 1 (right). Empty symbols correspond to non-breaking 
droplets, filled symbols to breaking droplets; diamonds correspond to simulations preformed with 
C = 1/15, while squares refer to simulations with C = 1/20. The solid lines show the separating 
boundary from equation Q29|) . with the viscosity dependent coefficients a(l/8) = 1 and a(l) = 0.5. 
Top: the frames in a) and b) show the motion of droplets with different elongation and same 
capillary number, corresponding to the points labeled a and b in the phase diagram for A = 1/8. 



binary mixture approaches that of two immiscible fluids only in the limit of an interface 
thickness which is vanishingly small compared to the droplet size. Also, we were able to 
show that the resistance of the droplet to breakup depends on the viscosity ratio of the two 
fluids, following the characteristic trend of isolated free droplets in a purely elongational 
flow. 

The main advantage of the method stems exactly from the fact that surface stresses are 
distributed over the finite interface thickness, allowing coarser meshes and therefore shorter 
simulation times. Furthermore, compared to other diffuse-interface models, the transport 
equations are derived from a generalized free-energy functional, which determines also the 
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equilibrium properties of the system, and is the starting point for extending the model to 
include more components and phases. The efficiency of the method, which we have shown 
does not come at the expense of accuracy or reliability, offers the possibility of handling 
two-phase flow in complex microfluidic devices; at the fundamental level, this should have 
an impact in the understanding of phenomena related to pattern formation in such devices, 
and we may expect that implementations of the method will prove to be useful for design 
purposes. 
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APPENDIX: SPURIOUS CURRENTS 



The origin of spurious currents can be explained by analyzing how the discretized form 
of differential operators modifies the continuum transport equations ()18|) and (j2*Uj) . Given 
set of stencil vectors {e^} and the lattice spacing Ax, the action of V a and V 2 on a field 
/(r) can be written in the general form 

V Q /(r) = — ^/(r + AzeOeta (A.l) 

i 

where the coefficients {ti} are chosen to ensure the equality between the discretized and 
continuum operators to the leading order. Their values are obtained by expanding in Taylor 
series the term /(r + Axe«) in equations (|A.1J) and (|A.2J) ; for a stencil having inversion 
symmetry, the expansion contains only the lattice tensors with an even number of indices 

— ^ "\ tjBja&ip, (A-3) 
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Tafias ~ ^ 1 ^ e ia e i/3 e «7 e i<5) (A-4) 



and the effect of the anisotropy of the lattice on the actual expression of the differential 
operators depends on the invariance properties of the tensors T^ n > . In the case of the 
simple cubic stencil with the 6 + 1 base vectors {e^} = (±1,0,0), (0, ±1,0), (0, ±1,0) and 
e = (0, 0, 0), the values t = for the gradient, t = —6 for the Laplacian, and ^=1,...^ = 1, 
imply that T^) = 25 a p is isotropic, such that V a /(r) = V Q /(r) + 0(Ax 2 ), and V 2 /(r) = 
V 2 /(r) ± 0(Ax 2 ). The coefficients of the fourth order tensor T^ 4 -*, on the other hand, 
depend on the orientation of the stencil vectors with respect to the coordinate axis, and are 
responsible for the anisotropy effects in the discretized transport equations. 

We now give a simple qualitative argument to relate the magnitude of the parasite currents 
to the dimensionless groups We, Re and C and the lattice spacing Ax. Let us consider 
the example of the equilibrated droplet in absence of any external flow. The transport 
equations (fT8|) and (}20|) are simplified to fi = const., V pP a p = and v = 0; in particular, 
the components F a of the volume force are 

F - = v - v + ¥m v °^ + ma^^ = (A - 5) 

or 

F a = V a p + -^VVV Q y? = 0, (A.6) 

where we have redefined the isotropic pressure in the interface region, p = p ± 
(| V<f\ 2 /2)(C/WeQ). The discretized form of the volumetric force F a is obtained with the 
replacement V a — ► V a , V 2 — > V 2 in equation (jA.6|) . When the fields and tp^ solving 
equation (|A.6|) for the equilibrated droplet are collocated on a grid, there is a residual stress 
of order 0(Ax 2 ) which depends explicitly on 



i^ (0 V 0) ] 

Ax 2 

±0(Ax 4 ) ^ (A.7) 



T (4) r (4) 

2-3! PlSP + WefT Q ^ } 4! Pj5eV 



and is balanced (neglecting the non-linear advection term in (|2UJl ) by the shear stress pro- 
duced by a non-zero velocity field v' 1 ', such that F a = V 2 v$ /Re. For the case of equal 
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viscosities, given that p(°' in the interface region is proportional to 1/ (WeC) and that each 
partial derivative introduces a factor 1/C, the module of the spurious currents is estimated 
by the simple relation 



and the proportionality factor depends on the non-invariant tensor which is responsible 
for the angular dependence of the parasite currents along the surface of the static droplet. 
Equation (jA.8|) includes the main factors governing the the magnitude of spurious currents, 
which are small for low surface tension or high viscosity, and can be controlled by including 
more grid points in the interface region, therefore decreasing the ratio Ax/C. This option is 
the most straightforward to use, at the cost of larger grids and therefore longer computational 
times. A more effective possibility consists in improving the symmetry properties of the 
discretization stencil {ej} by adding more vectors, such that the becomes isotropic and 
corrections are reduced to order (Ax/C) 4 
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